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ALGORITHMS AND MODELS FOR TURBULENCE NOT AT 
STATISTICAL EQUILIBRIUM 

NAN JIANG* AND WILLIAM LAYTON+ 


Abstract. Standard eddy viscosity models, while robust, cannot represent backscatter and have 
severe difficulties with complex turbulence not at statistical equilibrium. This report gives a new 
derivation of eddy viscosity models from an equation for the evolution of variance in a turbulent 
flow. The new derivation also shows how to correct eddy viscosity models. The report proves the 
corrected models preserve important features of the true Reynolds stresses. It gives algorithms for 
their discretization including a minimally invasive modular step to adapt an eddy viscosity code to 
the extended models. A numerical test is given with the usual and over diffusive Smagorinsky model. 
The correction (scaled by 10 —8 ) does successfully exhibit intermittent backscatter. 
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1. Introduction. Eddy viscosity models are the workhorses of practical turbu¬ 
lent flow simulations, |B\ Due to the wide experience with them, their limitations 
are also well recognized. They cannot represent backscatter (intermittent energy flow 
from turbulent fluctuations back to the mean velocity) without ad hoc fixes (called 
’’absurdities” in |23]) like negative viscosities. This report shows how to correct eddy 
viscosity models systematically to include backscatter based on a new and fundamen¬ 
tal derivation of eddy viscosity models. 

To begin, given an ensemble of initial conditions 

u(x, 0; ujj) = uo(x-, = 1, • ■ •, J, x £ fl, 
let u(x,t;ujj),p(x,t;ojj) be associated solutions to the Navier-Stokes equations (NSE) 

ut + u- S7u — vAu + Vp = f(x, £), and V ■ u = 0, in Q, (1-1) 

u = 0 on dfl. 


Let (•) denote ensemble averaging 

1 J 

( u ) ( x , t) := — '^2 u ( x > w i) and u\x, f; ujj) = u(x, t ; ujj) — (u ) (x, t). 

J 3 = 1 

Ensemble averaging the NSE yields the non-closed system: V • (u) = 0 and 

(u} t + (u) ■ V ( u) — V A (it) — V • R(u, it) + V ( p) = f{x, t), (1.2) 

where the Reynolds stress R(u, u) is 

R(u, it) := (it) ® (it) — (it (g> it) = — (i/ (g> u '), 

e.g., d, m, m■ Statistical models of turbulence begin with ensemble averaging and 
replace R(u, it) by an enhanced viscous term depending only on the mean velocity. 
We show in Section 2 that these eddy viscosity models are based on three steps. 
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1. The Boussinesq assumption (from [4j, [26], proven in [18]) that turbulent 
fluctuations (the action of V -R(u,u) in ( |1.2[ >) are dissipative on average in (1.2). This 
is followed by assuming that space and time averaged dissipativity holds pointuiise in 
time and space. 

2. The eddy viscosity hypothesis that this dissipativity aligns with the gradient 
or deformation tensor and thus can be represented by a viscous term with a turbulent 
viscosity coefficient v T ((u}), [251 . 

3. Model parametrization/calibration is done by fitting the turbulent viscosity 
coefficient vt({u)) to flow data. Calibration is equivalent to specifying a fluctuation 
model for Vu' in terms of V (u). 

The resulting eddy viscosity model (whose solution w(x,t), q(x,t ) is intended to 
be an approximation of the true flow averages (it), (p)) results: V • w = 0 and 


wt + w ■ S7w — V • ([i/ + ut{w)\ Vro) + Vq = f(x, t), in Q, (EV) 

w = 0 on dfl and w(x , 0) = (uq) in S2. 


Eddy viscosity models, with increasingly complex equations determining vt{w ), are 
the models of choice for most industrial turbulent flows, [B], and many parameteriza- 
tions of eddy viscosity models are known, e.g., [3], [5] [12], [13], [19], [20], [22], [31]. 
They have well recognized limitations in not modeling complex turbulence, backscat- 
ter or turbulence not at statistical equilibrium, e.g., 0, [21], [27], [30]. (The second 
assumption that the dissipativity of the Reynolds stress term aligns with V (u) also 
fails for some flows, [2T], but is not the issue addressed herein.) 

The correction required for eddy viscosity models to represent backscatter in 
non-statistically stationary turbulence, the case when the action of the fluctuations 
is intermittently non-dissipative, is derived and analyzed herein. Given the eddy 
viscosity parameterization vt(w), choose a re-scaling parameter f3 > 0 and define 

a(w) := \J 


The corrected EV model (derived in Section 2) is then V • w = 0 and 


Wt + (3 2 a(w )— (a(w)w) + w • Vu> 
at 

-V • ({v + v T {w)\ Vro) + Vq = f(x, t). 


(Corrected EV) 


In Section 3 time averaged dissipativity, an important feature of the true Reynolds 
stresses, is proven to be preserved in (Corrected EV). 

The (Corrected EV) differs from (EV) by the extra term /3 2 a(u>)^ ( a(w)w ). This 
term means time discretization introduces new issues, especially in adapting legacy 
codes from (EV) to (Corrected EV). Section 4 shows how time discretization can be 
done and preserve these important model properties, including the important case 
of modular adaptation of a legacy code available for (EV). Phenomenology is used 
to obtain some insight into calibration of the re-scaling parameter [3 in Section 5. 
Section 6 tests correction of the over-diffused Smagorinsky model. Even with a very 
small rescaling of the correction, /3 2 ~ O(10 -8 ), the numerical test shows that the 
corrected model does exhibit backscatter. 

2. Derivation of Corrected EV Models. Beginning with two results from 
[13] , this section shows that eddy viscosity models are based on three assumptions 
and that they are only consistent in a time averaged sense. Next we show how to 
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take a given eddy viscosity model and extend it to be energy consistent with the NSE 
pointwise in time. 

The ensemble averaged Navier Stokes equations, (1.2) above, involve the non- 
closed Reynolds stress R(u,u) = — (it'® it'). This term, which must be modelled, 
accounts for the effects of the fluctuations on the mean flow, e.g., 0 , i, m , ee 
[ 28]. Let || ■ ||, (•, •) denote the usual L 2 (Q) norm and inner product. Taking the inner 
product of (1.2) with (it) and performing the usual steps gives the equation for the 
kinetic energy evolution of (it): 

J t \\\( u ) II 2 + HI V ( u ) II 2 + X' R ( U ’ U ) : V = 

Thus, the effect of fluctuations on the mean flow is determined by the sign of 

RS(t) := I R(u,u) : V (it) dx. 

J n 

When RS(t) > 0, the effect of R(u,u) is dissipative while when RS(t) < 0, (volume 
averaged) backscatter occurs and fluctuations increase the energy in the mean flow. In 
EH], EH, two key properties of this Reynolds stress term were proven: time averaged 
dissipativity and an equation for the evolution of variance of fluctuations: 


LIM T ^l- [ RS(t)dt = [ [ v(\Vu'\ 2 )dxdt>0, (2.1) 

1 Jo 1 Jo Jn 

J R(u, u ) : V (it) dx = - — J ( u ' ■ u') dx + v J (Vi/ : \7u') dx. (2.2) 

Eddy viscosity models then follow from three assumptions. First, the statistical equi¬ 
librium assumption that dissipativity holds approximately at every instant in time 


f R(u,u) :V (u) dx ~ f v (|Vzt'| 2 ) dx. 
J n J n 


(2.3) 


The second is that Vu' aligns with V (u). Third, calibratiorQ provides a model of the 
fluctuations in terms of the mean flow 


action(Vu') ~ a((u))V (u). 


(2.4) 


Thus, 


/ R{u 7 u) : V (u) dx ~ / ^a((w)) 2 V (u) : V (it) dx 
Jn Jn 

= / —V • (t/a((n)) 2 V (it)) • vdx evaluated at v = (it). 

Jn 

Letting = va{{u )) 2 , this yields the eddy viscosity closure 

—V • R(u,u) <= — V • (t'j,((it))V (it)) + terms incorporated in Vp. 


Far from equilibrium, step (2.3) omits \ ^ J (u'j ■ it') dx in (2.2). This is the term 
that accounts for backscatter and other non-equilibrium effects. To model this term, 


1 Alternately, the dissipation in (2.3) is (t/|Vu , | 2 ) is replaced in the model by vr ( (u))\ V (u) | 2 . 
Ideally, then (IVu'l 2 ) ~ V ( u ) | 2 so a((u)) = ^/u~ 1 ux{{u)\ 
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v! must be expressed in terms of (u). For this, the simplest (explored herein) is to 
rescale (by /?, Section 4) the fluctuation model (2.4), yielding 


action(u') ~ (3a((u)) (u) 


This assumption yields 


J R(u , u ):V( u) dx 


1 d 

2 dt 


f3 2 a((u )) 2 1 (u) ( x,t)\ 2 dx + 


va{(u )) 2 |V (u) ( x,t)\ 2 dx , 


arising from an anisotropic time derivative /3 2 a((u))J^ (a (it)) (zz)) and an eddy vis¬ 
cosity term —V • (i / t((u))V (zz)) where ut{{u)) = ua((u)) 2 . This gives the closure 
model 


-V • R(u,u) ~ P 2 a({u))^ ( a((u )) (zz)) - V • (v T {u))V (zz)) 

and thus we have the corrected model: V • w = 0 and 
d 

w t + P 2 a{w)— (a(w)w) + w ■ Vw + Vp — V • {[v + zz T (u;)]Vu;) = / . (2.5) 

3. Analysis of The Corrected Smagorinsky Model. The classic example 
of an over-diffused model is the standard Smagorinsky model for which 

i^Tiw) = (C s S) 2 \Vw\, where C s ~ 0.1 ,6 = Ax. 

(Other numerical values of C s are also used, [25] Table 1.) Various fixes for it include 
van Driest damping (reducing near wall model dissipation) and Germano’s dynamic 
selection of /i = C 2 (x,t). These are often successful but the latter leads to negative 
values of /j that model backscatter but can induce instabilities. Often these are clipped 
(fi <= max{/i, 0}) eliminating backscatter being represented in the model. 

In this section we prove that the corrected Smagorinsky model preserves the prop¬ 
erty of the true turbulent fluctuations that on long time average they are dissipative, 
Theorem 3.2. For the Smagorinsky model we have 

vt(w) = ( C s S ) 2 |Viu|, and thus a(w) = v~ l ^ 2 C s 5\/\S7w\. 

Let (3 > 0, a(w) = z/ -1 / 2 C s c)^/|Vu;| and consider 

w t + j3 2 a(w)^ (a(w)w) + w ■ V w + Vp — V • {[v + jz t (w)]Vu;) =/, ) 

V • w = 0 in O x (0, oo), > (3.1) 

w = 0 on 90, w(x, 0) = w o(x), in O. J 

The model dissipation function for the effect of the model of the Reynolds stresses on 
the kinetic energy in the mean flow will be denoted by 

f d 

MD(t) := / (3 2 a(w)— (a(w)w) ■ w + v T {w)\S7w\ 2 dx. 

Jn dt 

The numerical tests in Section 6 establish (Figure 6.2) that, even for very small /3, 
the extra term does allow the sign of the model dissipation MD(t ) to fluctuate. We 
prove the time average of MD(t) is non-negative. 
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We assume 


G L°°(0,oo;L 2 (r!)) lW ; 0 (s) G L 2 (H), 
ty is a strong solution of (3.1). 


Lemma 3.1. For a(w) = ^ ^C s 8^/\S7w\ we have for all w G lT 1,3 (f2) 

||a('u>) u; || 2 < C'(fi)^ -1 ( C s 5 ) 2 11Vtu|| 3 3 . 


Proof. By Holder’s inequality and the Poincare-Friedrichs inequality 
\\a(w)w\\ 2 = v" 1 (C s 5) 2 I \Vw\\w\ 2 dx 

J o 

< v- 1 ( CJ) 2 ||«;|| 2 3||Vu;|| L 3 < C^p- 1 (C s S)' 2 ||V W || 3 3 . 


□ 

We prove next that MD(t) dissipates energy in the time average sense. 
Theorem 3.2. Let C s > 0,<5 > 0,/3 > 0 and w be a strong solution of \3.1ty . 
Then 


1 

T 



w , a(w)w G L°°(0, oo; L 2 (H)), 

\u + pt{w])\7w : S7wdx I dt < C < oo, C independent of T. 


(3.2) 

(3.3) 


The model also satisfies long term balance between energy input and dissipation: for 
any generalized limit LIM 

LIMt^ooj. J (/ [p + v T {w])\S7w\ 2 dx^dt = LIM T ^ 00 ^ J ( f,w)dt. 
Further 

1 f T 1 f T 

0< lim inf — / MD(t)dt< lim sup — / MD{t)dt< oo. 

T->oo T J g T -*oo T J q 


Proof. Take the inner product of the corrected Smagorinsky model with w. This 
yields 

\ < <lt [H W H 2 +^ 2 H a H u; ]] 2 ] + j n [v + VT(w)]\X7w\ 2 dx = (f 1 w). (3.4) 

Let 

F ■= ^;ll/llL“( 0 ,oo;L=(n)) and y(t) := ^ [\\w\\ 2 +/3 2 \\a(w)w\\ 2 ] . 

By Lemma 3.1, 

2 y(t) = ||u;|| 2 + /? 2 ||a(u>)u>|| 2 < C [v + vt{uj )] \S7w\ 2 dx 
< C ^||Viy|| 2 + (C s 8) 2 11 15 , 3 ) ■ 
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Since (f,w) < |||Vu>|| 2 + ^||/|| 2 i, y(t) thus satisfies 

y'(t) + ay(t) < F{< oo) for some a > 0. 


This implies y{t) £ L°°{ 0,oo) which is the first claimed a priori bound. This bound 
now implies that (f,w)(t) £ L°°(0,oo). Integrating (3.4) over [0,T] and dividing by 
T gives 


2 T 


+ f3 2 \\a(w)w\\ 2 ] (T) +-^ j ^ J [v + ^r(w])|Vw| 2 da:^ dt (3.5) 


= ff [lko|| 2 + /? 2 ||o(u;o)u;o|| 2 ] + 1 


(f,w)dt. 


Since ( f,w ) £ L°°( 0, oo) this implies (3.3) holds which implies as T —► oo limit 
inferiors and superiors exist. The a priori estimates (3.2), (3.3) imply that (3.5) takes 
the form 

0(1) + !^ ( y J n ^ + UT ^ w ^ Vw ^ dxS j dt = °^ + f J Q (f’ w ) dt - ( 3 - 6 ) 

Letting T — > oo implies long term balance between energy input and dissipation. 
Consider now the time average of MD(t). We have 




T 


0 Jfl 

t r«2 


d 

f) — ( a(w)w) ■ ( a(w)w ) + z't(w)|Vw| 


dxdt 


y^lK^HI 2 + ( c s $) 2 l|Vtu ||| 3 


dt 


= ff (IKM T )M T )II 2 - ||a(w(0))w;(0)|| 2 ) + 1 (C s S) 2 \\Vw\\ 3 L3 dt. 


From (3.5|, the right hand side equals 

r T 


f MD ^ dt= f J Q (f’ w ) dt ~ff o v\\\hv\\ 2 dt- ^ [||w|| 2 (T) - |K|| 2 ] , 


/ 0 

while from (3.6) we have 
r T 


fj MD ^ dt= f J Q (f’ w ) dt ~ f J q v\\Vw\\ 2 dt 

= 0 ( 1 ) - 0 ( 1 ) + 1 £ (C s 6) 2 MVHlia dt. 


Thus the limit inferior of the time average of MD(t ) exists and is non-negative. □ 

4. Time Discretization of Corrected Eddy Viscosity Models. This sec¬ 
tion presents three unconditionally stable, linearly implicit time discretizations of 
(Corrected EV). Method 1 and the modular Method 2 are first order accurate. Method 
3 is second order accurate. All three preserve the essential feature of time averaged 
dissipativity, Proposition 4.3. The second method shows how a legacy code for solving 
(EV) can be adapted to solve (Corrected EV). We suppress the spacial discretization 
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to reduce notation and focus on essential points. Let the timestep and associated 
quantities be denoted as usual by 


timestep = k, t n = nk, f n (x) := f{x,t n ), 
w n (x) = approximation to w(x, t n ), and a n :=a(w n ), utf := vt(w u ), n > 0, 

and set a -1 = a 0 = a(w°). We consider time discretization of (Corrected EV) under 
no-slip boundary conditions. The first method is: given w°, find w n satisfying 

,n„.,n+1 _ n— l„..n 


w n+1 — w ri 


2 


+ P z d 


a w 1 — a w 


k k 

+w n ■ Vw" +1 - V • (\v + i/£] Vw n+1 ) + V<f +1 = f n+1 (x) in fi, (Method 1) 
V • w n+1 = 0 in n, and w n+1 = 0 on dfl . 


This method is linearly implicit. We shall prove in Theorem 4.1 that it is uncondition¬ 
ally, nonlinearly stable. In linearly implicit methods (considered herein) nonlinearities 
are lagged to previous time levels (or extrapolated), so there is no difficulty if v T is 
determined by solving other nonlinear equations, such as in the k — e model, [5], |221 . 
Comparing the model, this method and the ensemble averaged NSE, we see that 

1. True effect of fluctuations on means: RS(t) = J Q R(u, u) : V ( u ) dx. 

2. Model: MD(t) = f Q (3 2 a(w)-^ (a(w)w) ■ w + VTi.w)\S7w\ 2 dx, 

3. Discrete version: MD n+1 = f Q . w n+i + t/ g|Vw n + 1 | 2 rfs. 

Theorem 4.1. (Method 1) is unconditionally, nonlinearly energy stable. For any 

N> 1 


-||t« JV || 2 + Vllo^-i^ll 2 

2 " " 2 11 11 


N-l 


+ (I \w n+1 - w n II 2 + P 2 \\a n w n+1 - a rl-1 u> n | | 2 ) 

n=0 2 

JV -1 . 

+k I [v + vJf] \S7w n+1 \ 2 dx 

n=0 Ja 

+ -/3 2 ||a-W|| 2 ) +fc£(r +1 ,t« n+1 ). 

' n =o 


,,0 112 


Proof. Multiply through by k, take the L 2 inner product of (Method 1) with 
w n+1 and use ( w n ■ V w n+1 , w n+1 ) = 0. This gives 


|K +1 || 2 - (■ w n+1 ,w n ) + 0 2 \\a n w n+1 \\ 2 - /3 2 (a n - 1 w n ,a n w n+1 ) 

+k f [v + vJf] \X7w n+1 \ 2 dx = k [ f n+1 ■ w n+1 dx. 

J J f2 

The second and fourth terms are treated by the polarization identity 

(w n+1 ,w n ) = i|K +1 || 2 + i|K|| 2 - \\\w n+1 - w n \\ 2 , 

(a n ~ 1 w n , a n w n+1 ) = ^\\a n w n+1 \\ 2 + l\\a n - 1 w n \\ 2 - ^\\a n w n+1 - a n - 1 w n \\ 2 . 

Collecting terms and summing from n = 0 to N — 1 finishes the proof. □ 

Since Theorem 4.1 is an energy equality we can identify various effects: 
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• Model kinetic energy = lllu^H 2 + \(5 2 \\a N 1 w Ar || 2 

• model dissipation = f n i/tf\Vw n+1 \ 2 dx 

• numerical diffusion = \ (||w; ra+1 — u; ra || 2 + fi 2 \\ a n w n+1 — a n ~ 1 w n || 2 ) . 
Rewriting the energy equality in an equivalent form gives 

^ (IK +1 || 2 - IKII 2 ) + ^IK +1 - w n \\ 2 + j/||Vu> n+1 || 2 + 

(^ (||a ra u> n+1 || 2 - ||a n_1 w n || 2 ) + ^-\\a n w n+1 - a^w 11 ]] 2 + f v%\\7w n+1 \ 2 dx} 
l 2A: 2k Jq j 

= (f n+ 1 ,w n+1 ). 

The first and last line arise from the terms in the usual backward Euler discretization 
of the NSE. The second line (bracketed) is an equivalent form of MD n+1 . 

Lemma 4.2. For (Method 1) we have 

MD n+1 = (||a n ui rl+1 || 2 - ||a" -1 u! ?l || 2 ) 

An, 

+ ^|ja n w n+1 -a n -V l || 2 + [ v?f\X/w n+l \ 2 dx. 

2fc Jn 

Dissipativity on time average holds for all three methods with the same manipu¬ 
lations of their discrete energy equality. We record it here for (Method 1). 

Proposition 4.3 (Time Averaged Dissipativity). Suppose sup 0<n<oo ||/(f n )|| < 
oo. Then, for T)v = NAt, 


1 


N 


lim inf — At V MD 

r N ->-oo T/v \ 

\ n=0 


,n+l 


> o. 


Proof. The proof is a discrete analog of the continuous case and will be omitted. 

□ 

4.1. Modular Correction of EV Models. Given a code that computes an 
approximation to the EV model 

w t +w- S7w — V • ([z/ + vt(w)\ Vw) + Vg = f(x, t). (4.1) 

Algorithm 4.2 presents a minimally intrusive, modular, postprocessor to solve: 


d 

w t + f3 2 a(w)— ( a(w)w ) + w ■ \7w — V • ([v + vt{w)\ S7w) + Vg = /( x, t). (4.2) 


Precise stability analysis requires a specific choice of the algorithm used to solve (4.1). 
For this we select the simple, linearly implicit, backward Euler method. 

Derivation and Consistency Error. To derive the modular postprocessor, 
rewrite (4.1) and (4.2) as, respectively, 


= ) and y\t) + /3 2 a{y)— (a(y)y) = f{t,y). 


The postprocessing given in Step 2 below suffices. 
Algorithm 4.4 (Method 2). Given y n ,y n ~ l , 





Step 1: Calculate Kemp by: Utcmp k = /(t n +i, Kemp)> 
Step 2 : Postprocess to obtain y n+1 from Kem^ by: 

[l + p 2 a(y n ) 2 ] y n+1 = Kemp + /3 2 a(KMK _1 )K . 


Eliminating y^n P from Step 2 shows that 2 / n+1 satisfies (with a n = a(y n )) 


temp 

,.n+l 


n n+ 1_ n-1 n 

■ A "——r y 

k 


= f(t n+1 ,y^ p ). 

Although this is close to Method 1, Kemp not y n+1 occurs in the RHS. The LHS is 
clearly a first order approximation to y'(t) + f3a(y)(a(y)y)'. The RHS, f(t n +i, K;mp)> 
is a first order approximation to f(t,y) provided Kemp ~ K +1 = O(k). Rearranging 
Step 2 gives 

,a{y n )y n+1 - a{y n ~ 1 )y n 


v n+1 
cf temp 


n+1 


= kip d a(y n )- 


= k \P 2 a{y) — (a(y)y)j + 0{k 2 ) = 0(k). 

Thus, Algorithm 4.4 is first order accurate approximation of y'(t) + (3 2 a(y)(a(y)y)' = 

Unconditional Stability of the Modular Algorithm. The utility of Algo¬ 
rithm 4.4 thus depends on its stability. This is now analyzed for its application to the 
corrected EV model. Algorithm 4.4 for (4.2) reads as follows. 

Algorithm 4.5. For n > 0, given w n 
Step 1 : Find w^J^p satisfying 


™?emp ~ wT ‘ 


Vm 


,n+1 


fc ' “■ temp 

-V • ([!/ + v$] V<£ p ) + = f n+ \x) in Cl, 

V • Kemp =0 inn, and Kemp = 0 ondCl . 

Step 2: Given w^^p, Qtemp find w n+1 , q n+1 satisfying 

[1 + p 2 (a n ) 2 ]w n+1 + Vq n+1 = Kemp + P 2 a n a n ~ 1 w n in ft 

V • w n+1 = 0 in Cl, and w n+1 = 0 on dCl . 


(4.3) 


We prove Algorithm 4.5 is unconditionally stable. 

Theorem 4.6. Algorithm 4-5 is unconditionally, nonlinearly energy stable. For 
any N > 1 


2 1 W 


,,N 112 


? 2 ||a iV -K JV || 2 


AT-l 

+ Eo (IK +1 - <mpl | 2 + IIKeLp - «K 2 ) 

n—0 

N-1 

n —0 




N-1 

+kY / (f n+1 


n—0 


W 


n+1 

temp 


)• 
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Proof. Take the L 2 inner product of Step 1 with vjff rr [ lp and follow the proof of 
Theorem 1. This gives 

\ I Kemp 11 2 - \ I K 11 2 + \ I Kemp - I ? ( Ste P 1 Energy) 

+k [ [u + ^}\Vw^ p \ 2 dx = k(f n+1 ,w^ P ). 

Jn 

Consider Step 2. Taking the L 2 inner product with w n+l gives 

|K +1 || 2 +/3 2 ||a”w n+1 || 2 = K + m >’ n+1 ) + ^(a n_1 w n ,aV +1 ). 


The two terms on the RHS are treated with the polarization identity 


(KempK +1 ) 


(a n ~ 1 w n ,a n w n+1 ) 


^Ketnpl|K^|K +1 || 2 -^IKrp-^ +1 || 2 , 

^\\a n w n+1 \\ 2 + ^\\a n - 1 w n \\ 2 - ^\\a n w n+1 - o T *- 1 «; n ||, 


giving 


^IK +1 || 2 + lp 2 \\a n w n+1 \\ 2 - ^\\a^w n \\ 2 + \\\w^ p - w n+1 \\ 2 

+ Pl\\a n w n + i_ a n -i w n \\ = I| Ke ^ p || 2 . 

Insert the LHS for |||Kempl| 2 i 11 (Step 1 Energy). This gives 

^IK + 1 II 2 - ^IKII 2 + \p 2 \\a n w n+1 \\ 2 - ^IK-KH 2 

+ 111^^-^ll 2 + llIKeK 1 ,, - ^ +1 N 2 + K-|| 2 

+k f \v + I/?] I Vw^/dx = k [ r +1 • v%£ p dx. 

J n Jo. 

The result now follows by summing over n. □ 

A Second Order Time Discretization. We present a second order method 
comprised of an IMEX combination of BDF2 and AB2 adapted to the new kinetic 
energy term. It shortens the notation considerably to denote the linear extrapolation 
of a variable <j) to t n+1 by <p* n+1 : 

(j )* n+ 1 := 2 cj) n — <j) n ~ 1 , for <j) = w,a, vt- 


The method is: given w°, w 1 , w 2 and w 3 (found by another method) find w n+1 for 
n > 3 satisfying 


3 w n+1 - 4 w n + 1 

2k 

j2 *n+i 3 a* n+1 w n+l - 4a* n w n + 

+/ a 2k 

+w* n+1 • Vw" +1 - V • {[v + 4” +1 ] Vw n+1 ) + Vg" +1 = f n+1 ( x) 
V • w n+1 = 0 in ft, and w n+1 = 0 on dQ . 


(Method 3) 


in 12, 
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Theorem 4.7. The method (Method 3) is unconditionally, nonlinearly, long time 
stable. For any N > 3 

j(lk"ll 2 + IK" + T) 

+ 7 (/3 2 ||o" v «i"l| 2 + 0 2 \\la’ N w N - a*~- 1 i« A '- I || 2 ) 

JV-1 1 

+ E {-P 2 \\a* n+1 w n+1 -2a* n w n +a* n - l w n ~ 1 \\ 2 

n —3 

+ h\w n+1 -2w n + w n ~ 1 \\ 2 + k f [p + v* T n+1 ] \Vw n+1 \ 2 dx} 

4 Jn J 

n-i i 

= E Kf n+1 ,w n+1 ) + 4(lk' 3 H 2 + l|2^ 3 - W 2 1| 2 ) 

n—3 

+ ^p 2 \\a* 3 w 3 \\ 2 + (3 2 \\2a* 3 w 3 - a* 2 u; 2 || 2 ). 

Proof. Take the L 2 inner product of (Method 3) with w n+1 and multiply through 
by k. This gives 

* (IK + 1 II 2 + K" +2 || 2 ) - * (|K || 2 + |K n+1 || 2 ) 

+ * (/3 2 ||a* ri+1 u; Il+1 || 2 + /3 2 ||2 a* n+1 w n+1 - a*V|| 2 ) 

~(/3 2 \\a* n w n \\ 2 + p 2 \\2a* n w n - 1| 2 ) 

+ ^/3 2 \\a* n+1 w n+1 - 2a* n w n + a *«- 1 w «- 1 || 2 

+ ]\\w n+1 -2u> n Pw^W 2 + k [ [u + v* T n+1 ] \Vw n+1 \ 2 dx = k(f n+1 ,w n+1 ). 

4 Jn 

Summing up above equality from n = 3 to n = N — 1 completes the proof. □ 

5. The Re-scaling Parameter /?. If vt{w) is well calibrated, then with a(w) 
= \Jvt(w)/v , we begin with the fluctuation model 

action(Vu') ~ a((tt))V (it). (5.1) 

A fluctuation model relating actionfu') to (u) is also needed. It is plausible but 
incorrect to begin with action(u') ~ a((u))(u). For reasons developed next, this 
relation must be rescaled (by a factor (3 « 1) yielding 

actionfu') ~ /3 ■ a({u)) ( u). (5-2) 


Definition 5.1. Let u denote an ensemble of realizations of the NSE with per¬ 
turbed initial data. The associated turbulent intensities of u and V u are, respectively, 


I{u) 


<IMI 2 > 

ll(«)ll 2 


and 7(Vm) 


(IIVu'H 2 ) 

11 V (u) || 2 ’ 
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The development of analytic insight into the scaling parameter 0 is based on 
phenomenology associating statistical means and fluctuations respectively with large 
and small spacial scales. Recall that, e.g., in, 0 , m , m, for fully developed 
turbulence, kinetic energy is concentrated in the large scales ((||u , || 2 ) << || ( u ) || 2 ) 
while energy dissipation is concentrated in the small scales ((HV-u'l] 2 ) >> ||V (it) || 2 ), 
0 ,0, [24] . On the scales where one is significant, the other is negligible. This implies 
that for fully developed turbulence 


I{u) « 1 « J(Vu). 


Beginning with the fluctuation models (5.1|, (5.2) we thus have 


NNN)NII 2 „ (INI 2 ) ,,, (INNl 2 ) „ IKN)v N II 2 

II Nil 2 "II Nil 2 IIVNII 2- IIVNII 2 

The quantities 

INN) Nil 2 . INN)v N II 2 
II Nil 2 IIVNII 2 


both represent (squares of) weighted averages of a((u)). If (as expected) these are of 
comparable magnitude, 0 << 1 (as expected) and we obtain the estimate 


-N) 

f(V«) 


« 1. 


(5.3) 


Predictions of Phenomenology. In 3d, for fully developed, homogeneous, 
isotropic turbulence an estimate of this quotient (and thus 0) can be calculated from 
the K41 theory and its predicted, time-averaged, energy density distribution E{k) ~ 
ct£ 2 / 3 k~ 5 / 3 , [2J]. Associate means with a length scale (typically the given mesh-width) 
<5. The well resolved scales and unresolved scales are then, respectively, n/L < k < ir/5 
and n/S < k < Tr/rj. rj = Re^ 3 / 4 L is the Kolmogorov micro-scale. We then calculate 


, / S -2/3 

i:/s v E(k)dk (f)_-d) 

s:!lmdk (f)- 2/s 

0 and S « L 

■s\ 2 /3 /n 2/3 


In the limit Re 


^N- 

°°, V ~ 
I(u) 


. -2/3 


-(f) 


-2/3 ' 


:ir 


i -(if 3 

Similarly, we calculate for S » r) 

/;/; k 2 E( k) dk 


z) 


+ H.O.T.S. 


J(Vu) ~ 


= . s ) 


fw/L k 2 E(k)clk W 


4/3 


H.O.T.s. 


In particular, dropping higher order terms, I(u) — ( 5/L ) 2 ^ 3 << 1 << J(Vu) ~ 
( 5/rj ) 4 ^ 3 . Thus, to leading order, in 3d 


0~ 



Re- 1 ' 2 
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Remark 5.2 (The 2d Case). Phenomenology of two dimensional turbulence is 
more complicated. The simplest case for forced turbulence is when the model incorpo¬ 
rates some extra mechanism to extract energy from the largest scales, energy is injected 
in an intermediate scale and 5/L is smaller than the injection scale but much larger 
than the micro-scale. Adapting the above spectral calculation to 2d gives 


Mesh Dependence. Apparently one option is to calculate the turbulent intensi¬ 
ties on a given mesh and use these to find the re-scaling parameter ft. Unfortunately, 
the result is severely limited by the chosen mesh as we now develop. Suppose spacial 
discretization is performed by a standard, conforming finite element method based 
on a mesh of elements e with element diameter (the local mesh width) denoted h e . 
For meshes satisfying an angle condition eliminating nearly degenerate elements and 
piecewise polynomial finite element velocities Uh, the inverse property 

||Vu h ||£ 2 (e) < CiNvhf 1 \\uh\\L 2 (e) (5-4) 

holds, where the constant depends only on local polynomial degree and mesh geometry 
(element angles). This implies that 

||Vit/j|| < C7Ary/i~ 1 ||it/,|| where h := min e h e . 


Let Uh denote an ensemble of discrete velocities computed on the same, fixed 
mesh as described above. Define the characteristic length-scale L of the ensemble 
mean velocity (expected but not guaranteed to be large) by 


L - 1 


IIVMI 
II Mil 


Proposition 5.3. 


Suppose (5-4) holds. 


Then 


satisfies 


Ph ■= 


I{Uh) 

ny Uh ) 


CjNV Y < Ph < CpFjr- 


Proof. We have, by rearranging, 

ii.u h ) = (IKII 2 ) (I|VM,)II 2 ) = (IKII 2 ) 2 
J(v« h ) <iiv u ;n 2 ) ii mii 2 (iiwj 2 ) ' 

From the inverse estimate and the Poincare-Friedrichs inequality (noting that C FF 
= 0{diameter of f2)) we have 

Cpl\\u' h \\ < ||V<|| < C INV h- x \\u' h \\. 
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Thus, 


i—2 (IKII 2 ) 

PF (l|v<|| 2 ) 


5: 1 < cJNVh 


(IKII 2 ) 

(liv<|| 2 )' 


Rearranging, it follows from the definition of L that, as claimed, 


cj 2 v h 2 L~ 2 < 


I{u h ) _ (IKII 2 ) r _ 2 


I(S7u h ) 


|V<" 2 


)■ 


< c 2 pf l~ 


□ 

Remark 5.4. Phenomenology \5.3 | ) suggests that the true value of /3 is small, 
f 3 « 1. On an under refined mesh, a reasonable default choice of (3 seems to be 

P(x) = h e (x ) 2 or ft = ( min e h e ) 2 . 


6. Numerical Explorations. The corrected model and its time discretizations 
give a closure that is dissipative on time average. The remaining question is whether 
the correction incorporates backscatter, i.e., whether MD(t) changes sign while being 
nonnegative on time average. To test the theory, we consider the Smagorinsky model 
(rather than models which perform better in practical calculations). This is because 
the Smagorinsky model is over-diffused and among models in use likely the one for 
which backscatter would be most difficult to introduce. Next, since it is believed that 
an inverse cascade and backscatter are much more common in the physics of 2d flow 
at high Re than for 3 d flows, we have selected a 2d test problem. (It is also one 
for which we have done a number of detailed simulation of the evolution of velocity 
ensembles in m, el m, m- While not directly relevant herein, this experience 
with velocity ensembles for this flow was useful in validation.) 

Test Problem: 2D Flow Between Offset Circles. Pick 


= {(x,y) : x 2 + y 2 < r 2 and (x - Ci ) 2 + {y - c 2 ) 2 > rf}, 
r i = 3,r 2 = 0.1, c = (ci,c 2 ) = (^,0), 
f{x, y, t) = (—4y(l - x 2 - y 2 ), 4x(l - x 2 - y 2 )) T , 


with no-slip boundary conditions on both circles. The flow (inspired by the extensive 
work on variants of Couette flow, 0), driven by a counterclockwise force (with / = 
0 at the outer circle), rotates about (0,0) and interacts with the immersed circle. 
This induces a von Karman vortex street which re-interacts with the immersed circle 
creating more complex structures. This flow also exhibits near wall turbulent streaks 
and a central polar vortex that pulsates. We discretize in space using the usual finite 
element method with Taylor-Hood elements, uni, using the code FreeFEM++, El 
and in time using (Method 1). For the Smagorinsky model to simplify notations we 
have previously used the full gradient in vt- In the implementation, we have used 
Sij, the stain rate or deformation tensor, 

vt = (0.lAx) 2 |^|. 


The choice C s = 0.1 is common though not universal, see Table 1 in [35 . Here 
151 = \/2SijSij. All the theory of the previous sections applies to choosing S instead 
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Fig. 6.1. Shown above is the mesh used for the flow between two offest circles. 


of Vt«. We take Ax to be the length of the shortest edge of all triangles. We also 
take 


a(-) 



/3 = 8 * icr 5 , v = icr 4 , At = o.oi, t = 10 . 


Here T is the simulation time. The numerical solutions are computed on an under¬ 
resolved, Delaunay-generated triangular mesh with 80 mesh points on the outer circle 
and 60 mesh points on the inner circle, providing 18,638 total degrees of freedom, 
refined near the inner circle (see Figure 6.1). For this mesh the shortest edge of all 
triangles is min e h e = 0.0110964 and the longest edge max e h e = 0.108046. 

We compute the following quantities. 


MD = / p 2 a(t n ) 


a(t n )w(t n+1 ) - a(t n ~ 1 )w(t n ) 


At 


+ is T \\Vw n+1 \\ 2 , 


TMD = / /3 2 a(t n ) 


a(t n )w(t n+1 ) - a(t n ~ 1 )w(t n ) 


At 


w(t n+1 )da 


■w(t n+1 )dx 


EVD = [ v^\S7w n+1 \ 2 dx, 

Jn 

VD = v\\S7w n+1 \\ 2 . 


Note that if jd = 0 (i.e., if we were solving the usual Smagorinsky model) we would 
have MD = EVD > 0. Observe in the first plot of Figure 6.2 that with (3 = 8 * 10~ 5 , 
MD(t) is on time average positive (consistent with theoretical predictions) but there 
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are times when MD becomes negative and indicates backscatter. Thus the corrections 
to the eddy viscosity models do have built into them the possibility of representing 
backscatter. Various other statistics are also plotted in Figures 6.2 including TMD 
which represents the effect of the new term, the eddy viscosity dissipation term EVD 
and the viscous dissipation term VD. 


x 10 






Fig. 6.2. v = 1/10000, At = 0.01, m = 80, n = 60, /? = 8 * 10“ 5 . 
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7. Conclusions. We have shown that eddy viscosity models can be quite eas¬ 
ily adapted to non-equilibrium turbulence and incorporate backscatter without using 
negative turbulent viscosities. The modified eddy viscosity model has been tested 
successfully for the Smagorinsky model, chosen because it is over dissipative. Some 
preliminary and formal calculations were given for the scaling parameter /?. It may 
also happen that for accuracy a different fluctuation model may be needed for the 
viscous dissipation term and the kinetic energy term. This is an open question. 
Strong solutions of the new models have been proven to share the property of the 
true Reynolds stresses of being dissipative on time average. We have also given three 
methods for time discretization preserving this property, including a modular correc¬ 
tion for an eddy viscosity code. There are many important open questions including 
accuracy tests, existence of weak solutions to the new models, model calibration and 
extension to and testing for better EV models. 
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